Simulation of Flux Lines with Columnar Pins: 
Bose Glass and Entangled Liquids 
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Using path integral Monte Carlo we simulate a 3D system of up to 1000 magnetic flux lines by mapping it 
onto a system of interacting bosons in (2+l)D. With increasing temperature we find a first order melting 
of fiux lines from an ordered solid to an entangled liquid signalled by a finite entropy jump and sharp 
discontinuities in the defect density and the structure factor S{G) at the first reciprocal lattice vector. 
In the presence of a small number of strong columnar pins, we find that the crystal is transformed into 
a Bose glass phase with patches of crystalline order nucleated around the trapped vortices but with no 
overall positional or orientational order. This glassy phase melts into a defected entangled liquid through a 
continuous transition. 
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The seminal paper of Abrikosov in 1957 ||l| showed, 
within a mean field calculation, that for fields below a 
critical field Hc2 (T) a type II superconductor allows mag- 
netic flux to penetrate in the form of a triangular array of 
vortices each carrying a unit flux quantum 0o ~ hc/2e. 
This picture is modified considerably in the case of high 
Tc cuprates which are quasi two-dimensional and have 
high transition temperatures The most dramatic 

consequence of the enhanced thermal fluctuations on the 
Abrikosov lattice is its melting into a vortex liquid by 
a first order transition as indicated in resistivity, neu- 
tron scattering, specific heat, SQUID magnetometry and 
local Hall probe measurements. Peak effect Q mea- 
surements on NbSe2 have also found a reentrant phase 
diagram in the H-T plane . Theoretical calculations 
have primarily been confined to density functional ap- 
proaches on the layered system of pancake vortices |^ 
vahd for BSCCO, lattice simulations of frustrated 3D 
XY-models B, Monte Carlo or Langevin dynamics of 
the vortices B, and determination of the melting line 
using Lindemann criterion |p^ . Simulations of vortices, 
treated as line or polymer-like objects, have been done 
using an analogy with interacting bosons Q by varia- 
tional and path integral techniques |^,|ll| . 

The effect of pins on the structure and dynamics of the 
vortex lattice adds another dimension to vortex behavior. 
The original work of Larkin and Ovchinnikov |]l2| pre- 
dicted that for d < 4, point pins destroy the long range 
order of the vortex lattice and break it up into small 
ordered coherence volumes. From the point of view of 
practical applications where the aim is to minimize the 
motion of vortices in order to reduce dissipation, it is 
found that columnar pins along the vortex direction are 
extremely effective in pinning vortices [l^ . There have 
been very few simulations attempting to study the phases 
and phase transitions in the presence of columnar pins. 
The ones based on 3D XY models assume an underlying 



numerical grid which affects the vortex structures for the 
small system sizes and high vortex densities that have 
been studied ||l|,|l. 

With this motivation, in this Letter we present the 
first simulations of the structure and thermodynamics of 
up to 1000 vortices with and without strong columnar 
pins using a continuum path integral quantum Monte 
Carlo (PIMC) simulation of the corresponding interact- 
ing bosons in (2-|-l)D. There is no bias of a superimposed 
lattice and no prior input of the nature of wave functions 
or the type of the transition. We extend previous simula- 
tions 1 1 1 1 to include columnar pins on considerably larger 
lattices. 

Our main findings are as follows- In the pure system, 
we see the reentrant nature of the melting line in agree- 
ment with the predictions of Nelson [|| . We find a tri- 
angular solid phase at low temperatures T which melts 
into an entangled liquid phase at higher T by a first order 
transition with a finite entropy jump. The defect density, 
and structure factor at the first reciprocal lattice vector, 
5(G), also show sharp discontinuities at the transition. 
In the solid phase, S{G) ^ 0{N) and drops sharply to 
0(1) in the liquid phase. In the presence of 10% colum- 
nar pins, the vortex lattice transforms into a Bose glass 
phase at low T. The Bose glass consists of patches of or- 
dered regions with only short range positional and orien- 
tational order. This phase is qualitatively different from 
the Bragg glass phase with point pins [|l^ . With increas- 
ing temperature, there is a transition from a Bose glass 
to an entangled 'defected liquid' which is considerably 
smeared. 

Consider a system of N flux lines in 3D characterized 
by a bending energy ei = {(^o/(47rA)}^ |l^ where A is the 
London penetration depth. The inter-vortex interaction 
is approximated as V{rij) — eoKQ{rij/X) where eo — 
0q/ (Stt^A^) is the scale of the inter- vortex potential, Kq is 
the modified Bessel function of the first kind and ri{z)'s 
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are the 2D position vectors of the flux Hnes in a plane 
perpendicular to their length at a height z. Ko{r) ~ 
— ln(r) as r ^ and is valid in the high vortex density 
limit, the regime studied in Ref. [0. In this paper, we 
fix the density on the lower part of the reentrant melting 
curve and therefore retain the full Ko{r) form of the 
potential. 

The partition function of the flux lines is mapped onto 
the world lines or Feynman paths of bosons described 
by the Hamiltonian 
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The bosonic path integrals are in two spatial dimensions 
and one imaginary time direction, which is equivalent 
to the thickness of the sample along the direction of the 
magnetic field. All lengths are measured in units of A and 
energies in units of eo- The dimensionless de Boer param- 
eter A = kBT/{2ei\^eoy^^ is a measure of quantum fluc- 
tuations (the relative strength of the kinetic energy vs the 
interaction energy) of the boson system, or equivalently 
a measure of temperature T in the vortex system. Once 
the mapping is established, we do a PIMC simulation fl^ ] 
of interac ting bosons at an inverse boson temperature 
/3 = 1 /Tft [Q . Periodic boundary conditions are used in 
all directions and permutations between bosons are in- 
cluded If Tfc < A^p, where p = 7VAV(Area) is the 
dimensionless density of the bosons, exchanges between 
the particles are possible and can generate entangled vor- 
tex lines. 

In Fig. |l|(a) we show the structure of 1000 vortices 
at low temperature A = 0.045 where the vortex system 
shows a beautiful triangular lattice with small excursions 
of the vortex lines about their equilibrium positions seen 
in Fig. |^(a). For A > 0.064 the vortices melt into an en- 
tangled liquid as seen from the actual configurations of 
the vortex lines in Fig. ||(b). Concomitantly, the super- 
fluid density ps of the corresponding boson system jumps 
from zero in the solid phase to a finite value in the liquid 
phase. The entanglement correlation length which is 
the typical length along the z-direction for vortex lines to 
cross, is on the order of the sample size in the solid and 
drops to 4ao, where ao is the planar lattice constant, 
in the liquid phase. 

The degree of positional order is quantified by the 
structure factor S(k = G) at the first reciprocal lat- 
tice vector of the triangular lattice. As seen in Fig. ||, 
with increasing temperature there is a sharp transition 
from a crystalline to a liquid phase at Am ~ 0.064. In 
a temperature-driven first order transition, the free en- 
ergy is continuous across the transition but the internal 
energy and entropy are discontinuous. In analogy in the 
mapped boson problem at T = 0, the total energy is con- 
tinuous but there are sharp jumps in the potential and 
kinetic energies between A = 0.062 and 0.064. The melt- 
ing transition is located by the intersection of the poly- 



nomial fits to the energy in the liquid and solid phases 
at Am = 0.064 in good agreement with Ref. ^. 

In the crystal S{k) shows Bragg peaks with S{G) ^ 
0{N), where N is the number of vortices, while above 
A„, S{k) is hquid-like with ^(G) - 0(1). By fitting 
the structure factor data to the form S{G) = 1 + {N — 
1) exp (— G^(m^)) for a finite system we have calculated 
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the Lindemann number cl = {u^) /clq = 0.24 at the 
transition just within the solid phase, in good agreement 
with other works j^], where (u^) is the mean squared dis- 
placement of the particles from their lattice positions . 

In addition to the positional correlations, the ori- 
entational order for six-fold symmetry is described by 
^6(ri) = (1/6) s'^^'J , where the sum over j runs 
over all the nearest neighbors of i in a triangulation 
picture, and 9ij is the angle the bond between par- 
ticles at Ti and Vj makes with an arbitrary reference 
axis. The orientational correlation function defined as 
9e{T) = (^6(7')*6(0)) also shows long range order in 
the clean system at low T. We find that the correlation 
lengths for positional and orientational order can reliably 
be extracted only from data on large system sizes ~ 1000 
vortices and we find that both types of order vanish at 

Am- 

A measure of imperfections in the triangular arrange- 
ment of the flux lines is the number of topologically- 
defected non six-fold coordinated sites in the Delaunay 
triangulation |21| of imaginary time slices of the system. 
This defect density is small at low temperatures in the 
clean system and shows a sharp jump across the melting 
transition (Fig. 
Effect of Columnar Pins: 

In the boson picture, randomly placed potential wells 
are correlated in the imaginary-time direction and cor- 
respond to columnar pins in the vortex representation. 
We model the columnar pin as an attractive delta func- 
tion which traps a vortex all along its length. In the 
boson picture the trapped vortex corresponds to a fixed 
classical particle. The vortex phase in the presence of 
columnar pins has been previously referred to as a Bose 
glass phase Our simulations provide the first de- 
tailed description of the structure and properties of this 
so-called Bose glass. In Fig. 0(b) we show the configu- 
ration of the 1000 vortex system at A = 0.045 with 10% 
of the vortices, those trapped by columnar pins, held 
fixed. The pins are found to seed or nucleate an ordered 
crystalline patch, however, since the pins are randomly 
located, the different patch orientations are not commen- 
surate and frustrates the translational and orientational 
order of the lattice. Columnar pinning dramatically re- 
duces S{G) from around 250 in the pure 1000 vortex 
system to about 6 in the presence of pins at A = 0.045 
in the Bose glass phase (Fig. |](a)). In addition, S{k) has 
a jagged peak, with shoulders and split peaks, which are 
characteristic of glasses. The Bose glass is rather differ- 
ent from a structural glass or a frozen liquid which has 
a peak height of iS'(G) « 1.5. At higher temperatures 
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the Bose glass melts into a defected liquid phase where 
the unpinned vortices are entangled (Fig. p(c)). In the 
parameter region investigated, we do not find any ev- 
idence for a disentangled liquid either with or without 
pins. The Bose glass phase that we find in these simula- 
tions is in a different limit from that studied by Nelson 
and Vinokur [ p2[ where the number of columnar pins is 
much larger than the number of vortices. 

The topological defects generated by the columnar pins 
in the vortex lattice destroy the positional and orienta- 
tional order within few lattice spacings (Fig. |^(b)) in con- 
trast to a Bragg glass in the presence of weak point 
pins which retains long range orientational order and al- 
gebraic positional order. 

Columnar pinning has a marked effect on the first order 
melting transition. The kinetic and potential energies 
vary smoothly through the transition. The sharp drop in 
S{G) in Fig. also disappears. The sharp jump in defect 
density is replaced by a much smoother increase through 
the transition as seen in Fig. ^ The local variations in 
the melting temperature caused by the quenched random 
potential below a certain critical dimension can cause a 
sharp first order transition to get rounded off in presence 
of quenched randomness . It is possible that at lower 
pin densities, the first order jumps in various quantities 
persist with reduced magnitudes. Further simulations 
are in progress to study the evolution of the sharp first 
order transition with increasing defect density. 

P. Sen would like to acknowledge B. Militzer, G. Bauer 
and E. Draeger for help with the code. We would also 
like to thank S. Banerjee, S. Bhattacharya, C. Dasgupta, 
T. Giamarchi, A. Grover, G. Menon, and A. Paramekanti 
for useful discussions. We acknowledge the support of the 
NCSA, UIUC for computational resources and software 
and the NSF grant no. DMR 98-02373. 
* e-mail: ntrivedi@tifr.res.in 



[lo; 

[11 
[12; 

[is: 

[14 

[is: 
[16: 



[17: 
[is: 



[19 

[20 

[21 

[22: 
[23: 



and D. Stroud, Phys. Rev. B 54, 1320 (1996); Otterlo et. 
al. Phys. Rev. Lett. 84, 2493 (2000). 
A. Houghton, R. A. Pelcovits, and A. Sudbo, Phys. Rev. 
B 40 6763 (1989); G. Blatter and H. Nordborg, Phys. 
Rev. B 54, 72 (1996). 

H. Nordborg and G. Blatter, Phys. Rev. Lett. 79, 1925 
(1997); Phys. Rev. B 58, 14556 (1998). 
A. I. Larkin and Yu. V. Ovchinnikov, J. Low Temp. Phys. 
34, 409 (1979). 

L. Civale et. al. Phys. Rev. Lett. 67, 648 (1991); Ch. 
Simon et. al. Nucl. Instr. and Meth. in Phys. Res. B 107, 
384 (1996); G. Pasquini et. al. Physica C 274, 165 (1997). 
K. H. Lee, D. Stroud and S. M. Girvin, Phys. Rev. B 48, 
1233 (1993). 

T. Giamarchi and P. Le Doussal, Phys. Rev. B 52 1242 
(1995). 

D. R. Nelson in The Los Alamos Symp. 1991, Applica- 
tions of High-temperature Superconductors, Ed: K. S. Be- 
dell et. al. (Addison- Wesley Pub. Co. 1992). 
D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). 
We use /3 = l/Tt = 20000 which is discretized into M = 
128 slices. The time step 5t = (3/M = 156.25 and the 
typical energy E ~ 10~^, which ensures that the Trotter 
error 5tE <^ 1. 

It is possible to have an entangled liquid with either 
bosons or boltzmannons, but bosons generate the Tt = 
state at a higher temperature and are thus computation- 
ally more convenient [ p^ . 

The mean squared displacement (m^) can also have some 
finite size effects, see E. Draeger and D. M. Ceperley, 
Phys. Rev. B 61 12094 (2000). 

F. F. Preparata and M. L. Shamos, Computational Ge- 
ometry: An Introduction, (Springer- Verlag, New York, 
1985). 

D. R. Nelson and V. M. Vinokur, Phys. Rev. B 48 13060 
(1993). 

Y. Imry and M. Wortis, Phys. Rev. B 19 3580 (1979). 
FIGURE CAPTIONS 



[1] 
[2] 

[3] 



[4] 
[5] 

[6] 

[7] 
[8] 



[9] 



A. A. Abrikosov, Sov. Phys. JETP 5, 1174 (1957). 

G. Blatter et. al.. Rev. Mod. Phys. 66, 1125 (1994). 

H. Safar et. al. Phys. Rev. Lett. 69, 824 (1992); R. Cubitt 
et. al. Nature 365, 407 (1993); U. Welp et. al. Phys. Rev. 
Lett. 76, 4809 (1996); E. Zeldov et. al. Nature 375, 373 
(1995); A. Schilling et. al. Nature 382, 791 (1996); 

K. Ghosh et. al. Phys. Rev. Lett. 76 4600 (1996). 

D. R. Nelson, Phys. Rev. Lett. 60, 1973 (1988); D. R. 

Nelson and H. S. Seung, Phys. Rev. B 39, 9153 (1989). 

W. R. Magro and D. M. Ceperley, Phys. Rev. B 48, 411 

(1993). 

S. Sengupta et. al. Phys. Rev. Lett. 67, 1023 (1991). 

(a) R. E. Hetzel et. al. Phys. Rev. Lett. 69, 518 (1992); 

(b) Y. Li and S. Teitel, Phys. Rev. B 47, 359 (1993) and 
49, 4136 (1994); (c) A. E. Koshelev, Phys. Rev. B 56, 11 
201 (1997); (d) S. Ryu and D. Stroud, Phys. Rev. Lett. 
78, 4629 (1997); (e) X. Hu, S. Miyashita and M. Tachiki 
iMd. 79, 3498 (1997). 

S. Ryu et. al. Phys. Rev. Lett. 68, 710 (1992); S. Ryu 



FIG. 1. Time averaged locations of 1000 vortices at 
A = 0.045 on a slice perpendicular to the z-direction and 
then averaged over many slices along the thickness of the vor- 
tex system; red regions represent higher vortex density and 
blue regions are lower vortex density, (a) Triangular lattice of 
vortices in the pure system, (b) Frozen Bose glass phase with 
10% strong columnar pins for one realization of disorder. A 
vortex, with a delta function density distribution, is trapped 
at the center of a pin (not shown). The blue region is the re- 
gion of reduced vortex density from which other vortices are 
repelled by the trapped vortex. 
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FIG. 2. A snapshot of vortex positions on the top layer 
(shown as black dots) and the projections on to the top layer 
of the vortex line (shown in red). Only a portion of the 1000 
vortex system is shown for clarity. (a)For A = 0.045 in the 
solid phase there arc very small excursions of the vortex lines. 
(b)For A = 0.065 just above the melting transition in the clean 
system the vortex lines are in a highly entangled liquid phase 
and the vortex density upon averaging is uniform. (c)For 
A — 0.065 with 10% pins the system is also melted into a 
defected entangled liquid. In the portion shown there are 3 
pins each trapping one vortex, (isolated black dots), which 
excludes other vortices around it. 

FIG. 3. The structure factor ^(k = G) at the first recip- 
rocal lattice vector of the triangular lattice as a function of 
A for 100 vortices at a density p = 0.02. The circles show 
S{G) in the clean system which shows a transition from a 
crystalline to a liquid phase at A ~ 0.064. The squares show 
^(G) with 10% columnar pins which shows a much reduced 
degree of translational order in the Bose glass phase. 

FIG. 4. Fraction of non six-fold coordinated sites in the 
clean and pinned systems. In the clean system the defect 
density shows a sharp jump across the transition at A = 0.064. 
In the presence of pins, the topological defect density is higher 
within the phases and the transition is smoothed out. 

FIG. 5. A 1000 vortex system with 10% columnar pins at 
A = 0.045 in the Bosc glass phase, (a) The structure factor 
S{k) vs. k showing a much reduced peak from about 250 in 
the crystal to ~ 6 in the Bose glass phase; but higher than 
in the liquid with peak height 1.5. (b) The radial positional 
distribution function g(r) and the orientational correlation 
function ge{r). An envelope of the form exp(r/^) can be fit 
for both g{r) and geir) with ^ = 13A and f = 15 A respectively. 
For density p = 0.02, the inter-vortex separation ao = 7.6A so 
the decay length of the correlations is ~ 2ao in the Bose glass 
phase. 
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This figure "figl.gif" is available in "gif" format from: 



http://arXiv.org/ps/cond-mat/0008405vl 
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